#############################################################
#
# File:     Maps.R
# Project:  Partisan Polarization on Black Suffrage, 1785-1868
# Function: Creates maps of Free Black Population by County
#           and of various measures of anti-slavery support
#           (i.e., Figures 2 and 5)
# Author:   David A. Bateman
# Date:     February 26, 2019 (this version)
#
#############################################################
rm(list=ls(all=TRUE))
R.Version()

library(shapefiles)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(ggplot2)
library(dismo)

dir <-paste("d:/dab/Dropbox/Documents/CSDP/Black Suffrage/", sep="")
setwd(dir)


#####################################################################################################
##
## Figure 2

## Input the shape files and population data
input1 <- paste("FinalFiles/Shapefiles/Counties1840.shp", sep="")
input2 <- paste("FinalFiles/data/FPOC1840.dta", sep="")
input3 <- paste("FinalFiles/Shapefiles/States1840.shp", sep="")
## Input population data for MidAtlantic inset.
input4 <- paste("FinalFiles/data/FPOC1840Midatlantic.dta", sep="")

states.map <- readShapeSpatial(input3, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
fpoc <- read.dta(input2)
fpocMA <- read.dta(input4)
names(fpoc)[1]<-"FIPS"
names(fpocMA)[1]<-"FIPS"

## Clean population files
fpoc <-fpoc[,c("FIPS", "black1843"),drop=FALSE]
fpocMA <-fpocMA[,c("FIPS", "black1843"),drop=FALSE]


## Drop Unorganized Territory from Map
states.map <-states.map[states.map$TERR_TYPE!="Unorganized Territory",]

## Put the data into the polycounties file
polycounties <- readShapePoly(input1, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
polycounties <- SpatialPolygonsDataFrame(polycounties , data=as(polycounties , "data.frame"))

## Some of the shapefiles have the wrong FIP codes. Fix these here 
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45075")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45013")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45043")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45019")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45069")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45051")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45089")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45029")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45075")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45076")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45011")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45085")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"45085")
levels(polycounties$FIPS)<-c(levels(polycounties$FIPS),"11000")
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==2578]<-11000 # DC
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==13085]<-45075
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==12879]<-45013
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==12987]<-45043 # Georgetown
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==12909]<-45019 # Charleston
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==12861]<-45051 # Horry
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==13108]<-45089 # Williamsburg
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==13088]<-45075 # Orangeburg
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==13094]<-45076 # Pendleton
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==12872]<-45011 # Barnwell
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==13177]<-45085 # Sumter
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==12910]<-45069
polycounties$FIPS[is.na(polycounties$FIPS) & polycounties$ID_NUM==12944]<-45029 # Colleton
polycounties$FIPS[polycounties$ID_NUM==12921]<-45057

### Set random distribution at county level
polycounties@data$order <- 1:nrow(polycounties@data)
polycounties@data <- merge(polycounties@data,fpoc, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
polycounties@data<-polycounties@data[order(polycounties@data$order),]

polycounties@data$order <- 1:nrow(polycounties@data)
polycounties@data <- merge(polycounties@data,fpocMA , by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
polycounties@data<-polycounties@data[order(polycounties@data$order),]

plotvar1 <- polycounties@data$black1843.x/10
plotvar1[is.na(plotvar1)] <- 0
dots.rand1 <- dotsInPolys(polycounties, as.integer(plotvar1), f="random")

plotvar2 <- polycounties@data$black1843.y/10
plotvar2[is.na(plotvar2)] <- 0
dots.rand2 <- dotsInPolys(polycounties, as.integer(plotvar2), f="random")

## Identify disenfranchising states
disfr<-c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1)
states.map@data$disfr <- disfr
states.map@data$color <- "grey40"
states.map@data$color[states.map@data$disf == 1]<-"white"
color.palette =  colorRampPalette(c("grey40","white"),space="rgb", bias=1, interpolate="spline")
cols <-color.palette(2)
index = c(cols[1], cols[2])

################################################################################
## Create Figure 2
label<-paste("FinalFiles/Figures/Figure2.pdf", sep="")
#label<-paste("FinalFiles/Figures/Figure2.eps", sep="")
#cairo_pdf(filename =label, width=8, height=6, onefile=TRUE, bg="transparent", family="serif")
pdf(label, width=8, height=6, bg="transparent", family="serif" )
par(mfrow=c(1,1),oma=c(1,0,1,0), mar=c(0,0,2,0))

plot(states.map, col=states.map@data$color)
plot(dots.rand1, add=T, pch=20, cex=0.025, col="black")

title(main = "Figure 2: Geographic Distribution of Free Black Americans \n and Disenfranchising States, 1843")
legend('left', bty="n", fill=index , density=c(NA, NA), ncol=1,
       x.intersp = 1, xjust=0, yjust=0, text.width=10,
       legend = c(	"States Recognizing Voting \nRights of Free Men of Color\n", "States Disenfranchising \nFree Men of Color"))

subtitlabel<-paste("Each dot is the approximate county location of ten free persons of color. \nShapefiles from Siczewicz (2011) and data from Haines and ICPSR (2010).")
mtext(subtitlabel, outer=TRUE, side=1, cex=1, adj=0, family="serif")

### Create inset on mid-Atlantic
int <-crop(states.map, dots.rand2)
shapefile(int, 'Analyses/Shapefiles/MidAtlantic.shp', overwrite=TRUE)

par(fig = c(0.65, .97, 0.03, 0.43), new = T)
plot(dots.rand2,  pch=20, cex=0.005, col="black")
plot(int, col=int@data$color, add=TRUE)
plot(dots.rand2,  pch=20, cex=0.005, col="black", add=TRUE)
box(lty = 'solid')
#dev.copy(postscript, label, horizontal = FALSE)

dev.off()

#############################################################################################
#############################################################################################
### Figure 5
#############################################################################################
#############################################################################################

### Load shapefiles and data. StatesXXX has states from years when votes were held
### input1 will be redefined with each panel of Figure
input2 <- paste("FinalFiles/data/MapsVotingData.dta", sep="")
input3 <- paste("FinalFiles/ShapeFiles/StatesXXXX.shp", sep="")

bsuff <- read.dta(input2)
state.map <- readShapeSpatial(input3, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))

label<-paste("FinalFiles/Figures/Figure5.pdf", sep="")
#cairo_pdf(filename =label, width=11, height=8.5, onefile=TRUE, bg="transparent", family="serif")
pdf(label, width=11, height=8.5, bg="transparent", family="serif")
par(mfrow=c(2,2),mar=c(0,0,2,0), oma=c(1.5,0,2,0))

#########################################################
### Create a function to draw heat maps
plot.heat1 <- function(counties.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  ##Break down the value variable
  if (is.null(breaks)) {
    breaks=
      seq(
        floor(min(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,
        ceiling(max(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,.1)
  }
  counties.map@data$zCat <- cut(counties.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(counties.map@data$zCat)
  if (is.null(col.vec)) col.vec <- color.palette1(length(levels(counties.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(counties.map@data$zCat) <- cutpointsColors
  plot(counties.map,border=gray(.9), lwd=bw,axes = FALSE, las = 1,col=as.character(counties.map@data$zCat), add=TRUE)
  
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
}
####################################################################################
##
## Figure 5(a)
input1 <- paste("FinalFiles/ShapeFiles/SuffrageCounties.shp", sep="")
counties.map <- readShapeSpatial(input1, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))

names(bsuff)[1]<-"FIPS"

voteshare<-bsuff[,c("FIPS", "bsuff"),drop=FALSE]
## get the order of the original counties.map@data
counties.map@data$order <- 1:nrow(counties.map@data)
counties.map@data <- merge(counties.map@data,voteshare, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
counties.map@data<-counties.map@data[order(counties.map@data$order),]
names(counties.map@data)[ names(counties.map@data)=="bsuff"]<-"noise"

color.palette1 =  colorRampPalette(c("white","black"),space="rgb", bias=1, interpolate="spline")
cols <-color.palette1(100)
test = c(cols[1], cols[51],cols[100], "black")

plot(state.map, density=c(30))
plot(counties.map, add=TRUE)
plot.heat1(counties.map,z="noise",breaks=c(0:100), plot.legend=FALSE)
plot(state.map, add=TRUE)
title(main = "Legislator Support for Black Suffrage, 1830-1860")
legend('bottomleft', bty="n", fill=test, density=c(NA, NA, NA, 30), ncol=2,
       x.intersp = 1, xjust=0, yjust=0, text.width=10,
       legend = c(	"Consistent\nOpposition", "Inconsistent or\nPartial Support",
                   "Consistent\nSupport", "No Votes"))
####################################################################################
##
## Figure 5(b)
input1 <- paste("FinalFiles/ShapeFiles/Liberty1844.shp", sep="")
counties.map <- readShapeSpatial(input1, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))

names(bsuff)[1]<-"FIPS"

voteshare<-bsuff[,c("FIPS", "liberty"),drop=FALSE]
## get the order of the original counties.map@data
counties.map@data$order <- 1:nrow(counties.map@data)
counties.map@data <- merge(counties.map@data,voteshare, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
counties.map@data<-counties.map@data[order(counties.map@data$order),]
names(counties.map@data)[ names(counties.map@data)=="liberty"]<-"noise"


color.palette1 =  colorRampPalette(c("white","black"),space="rgb", bias=1, interpolate="spline")
cols <-color.palette1(100)
test = c(cols[1], cols[25], cols[50], cols[75], "black")

plot(state.map, density=c(30))
plot(counties.map, add=TRUE)
plot.heat1(counties.map,z="noise",breaks=c(0,2.5,5,10,100), plot.legend=FALSE)
plot(state.map,add=TRUE)
title(main = "Liberty Vote Share, 1844")
legend('bottomleft', bty="n", fill=test, density=c(NA, NA, NA, NA, 30), ncol=2,
       x.intersp = 1, xjust=0, yjust=0, text.width=10,
       legend = c(	"0 - 2.5%", "2.5 - 5%",
                   "5 - 10%", ">10%", "No Votes"))

####################################################################################
##
## Figure 5(c)
input1 <- paste("FinalFiles/ShapeFiles/Freesoil1848.shp", sep="")
counties.map <- readShapeSpatial(input1, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))

names(bsuff)[1]<-"FIPS"

voteshare<-bsuff[,c("FIPS", "freesoil"),drop=FALSE]
## get the order of the original counties.map@data
counties.map@data$order <- 1:nrow(counties.map@data)
counties.map@data <- merge(counties.map@data,voteshare, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
counties.map@data<-counties.map@data[order(counties.map@data$order),]
names(counties.map@data)[ names(counties.map@data)=="freesoil"]<-"noise"


color.palette1 =  colorRampPalette(c("white","black"),space="rgb", bias=1, interpolate="spline")
cols <-color.palette1(100)
test = c(cols[1], cols[17], cols[35], cols[51], cols[85], "black")

plot(state.map, density=c(30))
plot(counties.map, add=TRUE)

plot.heat1(counties.map,z="noise",breaks=c(0,5,15,25,50,100), plot.legend=FALSE)
plot(state.map, add=TRUE)
title(main = "Free Soil Vote Share, 1848")
legend('bottomleft', bty="n", fill=test, density=c(NA, NA, NA, NA, NA, 30), ncol=2,
       x.intersp = 1, xjust=0, yjust=0, text.width=10,
       legend = c(	"0 - 5%", "5 - 15%", "15 - 25%",
                   "25 - 50%", "50 - 100%", "No Votes"))
###############################################################
##
## Figure 5(d)
input1 <- paste("FinalFiles/ShapeFiles/Republican1856.shp", sep="")
counties.map <- readShapeSpatial(input1, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))


voteshare<-bsuff[,c("FIPS", "republican"),drop=FALSE]
## get the order of the original counties.map@data
counties.map@data$order <- 1:nrow(counties.map@data)
counties.map@data <- merge(counties.map@data,voteshare, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
counties.map@data<-counties.map@data[order(counties.map@data$order),]
names(counties.map@data)[ names(counties.map@data)=="republican"]<-"noise"

color.palette1 =  colorRampPalette(c("white","black"),space="rgb", bias=1, interpolate="spline")
cols <-color.palette1(100)
test = c(cols[1], cols[35],  cols[65], cols[100], "black")

plot(state.map, density=c(30))
plot(counties.map, add=TRUE)

plot.heat1(counties.map,z="noise",breaks=c(0,40,50,65,100), plot.legend=FALSE)
plot(state.map, add=TRUE)
title(main = "Republican Vote Share, 1856")
legend('bottomleft', bty="n", fill=test, density=c(NA, NA, NA, NA, 30), ncol=2,
       x.intersp = 1, xjust=0, yjust=0, text.width=10,
       legend = c(	"0 - 40%", "40 - 50%",
                   "50 - 65%", "65 - 100%", "No Votes"))

mtext("Figure 5: Geographic Distribution of Black Suffrage Support and Antislavery", outer = TRUE, cex = 1)
mtext(side=1, "Shapefiles from Siczewicz (2011) and public voting data from ICPSR (1999).", outer=TRUE)

dev.off()